/*==============================================================================
FILE NAME: Figure_A4.do
CREATED: 13 June 2025
==============================================================================*/

**Figure A4

/* Set directory if working independently through code
if c(username)=="" { //insert username
	global rootdir "" // insert root path
	global processed_data "$rootdir/processed_data" 
	global figures "$rootdir/output/figures"  // Define global paths for replication package
} 
*/

// Panel A
// Load in data
use "$processed_data/Air_panel_with_referred.dta", clear

// Drop observations wtih no air investigations
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
drop if p_air_inv_`h' == .
}

// Generate placebo variables for months prior to investigation
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
drop if p_air_inv_neg`h' == .
}

// Initialize storage variables for plotting
cap drop b u d se b_ref se_ref u_ref d_ref Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0
gen b_ref = 0
gen se_ref = 0
gen u_ref = 0 
gen d_ref = 0

// Run regressions for 0 to 12 months after incident
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_inv_`h' p_referred p_air_incident, absorb(RN_year t) cluster(RN_id )
unique t if e(sample)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
replace b_ref = _b[p_referred] if Years == `h'
replace se_ref = _se[p_referred] if Years == `h'
replace u_ref = (_b[p_referred] + 1.96*_se[p_referred]) if Years == `h'
replace d_ref = (_b[p_referred] - 1.96*_se[p_referred]) if Years == `h'
}
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_inv_neg`h' p_referred p_air_incident, absorb(RN_year t) cluster(RN_id)
unique t if e(sample)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'
replace b_ref = _b[p_referred] if Years == -`h'
replace se_ref = _se[p_referred] if Years == -`h'
replace u_ref = (_b[p_referred] + 1.96*_se[p_referred]) if Years == -`h'
replace d_ref = (_b[p_referred] - 1.96*_se[p_referred]) if Years == -`h'
}
preserve
keep if Years != .
keep b u d se b_ref se_ref u_ref d_ref Years Zero

restore

bys RN_id: egen ever_referred = max(p_referred)
keep if ever_referred == 1

cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0


// Regressions for post-incident outcomes
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_inv_`h' p_air_incident, absorb(RN_year t) cluster(RN_id )
unique t if e(sample)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
}

// Regressions for pre-incident placebo outcomes
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_inv_neg`h' p_air_incident, absorb(RN_year t) cluster(RN_id)
unique t if e(sample)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

}


preserve
keep if Years != .
keep b u d se Years Zero
graph set window fontface "Times New Roman"
twoway(rarea u d Years, col(navy) fint(inten20) lwidth(0) lpattern(solid))(line b Years, lcolor(navy) lpattern(solid) lwidth(medium))(line Zero Years, lcolor(black)), xlabel(-12(1)12, nogrid labsize(vlarge)) ylabel(0(0.1)0.6, labsize(vlarge)) legend(off) ytitle("{&Delta} P(Investigation)", size(vlarge)) xtitle("Month", size(vlarge)) graphregion(color(white)) plotregion(color(white)) xsize(8.6)

// Save figure
graph export "$figures/Figure_A4_Panel_A.pdf", replace
restore

// Panel B
// Import data
use "$processed_data/Air_Panel_with_referred.dta", clear

//Drop missings
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
drop if p_air_inv_`h' == .
}
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
drop if p_air_inv_neg`h' == .
}

cap drop b u d se b_ref se_ref u_ref d_ref Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0
gen b_ref = 0
gen se_ref = 0
gen u_ref = 0 
gen d_ref = 0

// Regressions for post-incident outcomes
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_`h' p_referred p_air_incident, absorb(RN_year t) cluster(RN_id )
unique t if e(sample)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
replace b_ref = _b[p_referred] if Years == `h'
replace se_ref = _se[p_referred] if Years == `h'
replace u_ref = (_b[p_referred] + 1.96*_se[p_referred]) if Years == `h'
replace d_ref = (_b[p_referred] - 1.96*_se[p_referred]) if Years == `h'
}

// Regressions for pre-incident placebo
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_neg`h' p_referred p_air_incident, absorb(RN_year t) cluster(RN_id)
unique t if e(sample)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'
replace b_ref = _b[p_referred] if Years == -`h'
replace se_ref = _se[p_referred] if Years == -`h'
replace u_ref = (_b[p_referred] + 1.96*_se[p_referred]) if Years == -`h'
replace d_ref = (_b[p_referred] - 1.96*_se[p_referred]) if Years == -`h'
}
preserve
keep if Years != .
keep b u d se b_ref se_ref u_ref d_ref Years Zero

restore


bys RN_id: egen ever_referred = max(p_referred)
keep if ever_referred == 1

cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

// Regressions for post-incident outcomes
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_`h' p_air_incident, absorb(RN_year t) cluster(RN_id )
unique t if e(sample)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
}

// Regressions for pre-incident outcomes 
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_nov_neg`h' p_air_incident, absorb(RN_year t) cluster(RN_id)
unique t if e(sample)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

}


preserve
keep if Years != .
keep b u d se Years Zero
graph set window fontface "Times New Roman"
twoway(rarea u d Years, col(navy) fint(inten20) lwidth(0) lpattern(solid))(line b Years, lcolor(navy) lpattern(solid) lwidth(medium))(line Zero Years, lcolor(black)), xlabel(-12(1)12, nogrid labsize(vlarge)) ylabel(0(0.02)0.08, labsize(vlarge)) legend(off) ytitle("{&Delta} P(Notice of Violation)", size(vlarge)) xtitle("Month", size(vlarge)) graphregion(color(white)) plotregion(color(white)) xsize(8.6)

// Save figure
graph export "$figures/Figure_A4_Panel_B.pdf", replace
restore

// Panel C
// Import data
use "$processed_data/Air_Panel_with_referred.dta", clear

// Drop missings
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
drop if p_air_inv_`h' == .
}
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
drop if p_air_inv_neg`h' == .
}

cap drop b u d se b_ref se_ref u_ref d_ref Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0
gen b_ref = 0
gen se_ref = 0
gen u_ref = 0 
gen d_ref = 0

// Regressions for post-incident outcomes
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_noe_`h' p_referred p_air_incident, absorb(RN_year t) cluster(RN_id )
unique t if e(sample)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
replace b_ref = _b[p_referred] if Years == `h'
replace se_ref = _se[p_referred] if Years == `h'
replace u_ref = (_b[p_referred] + 1.96*_se[p_referred]) if Years == `h'
replace d_ref = (_b[p_referred] - 1.96*_se[p_referred]) if Years == `h'
}

// Regressions for pre-incident placebo outcomes
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_noe_neg`h' p_referred p_air_incident, absorb(RN_year t) cluster(RN_id)
unique t if e(sample)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'
replace b_ref = _b[p_referred] if Years == -`h'
replace se_ref = _se[p_referred] if Years == -`h'
replace u_ref = (_b[p_referred] + 1.96*_se[p_referred]) if Years == -`h'
replace d_ref = (_b[p_referred] - 1.96*_se[p_referred]) if Years == -`h'
}
preserve
keep if Years != .
keep b u d se b_ref se_ref u_ref d_ref Years Zero

restore


bys RN_id: egen ever_referred = max(p_referred)
keep if ever_referred == 1

cap drop b u d se Years Zero
gen Years = _n-13 if _n<=25
gen Zero = 0 if _n <=25
gen b = 0
gen se = 0
gen u = 0 
gen d = 0

// Regressions for post-incident outcomes
foreach h in 0 1 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_noe_`h' p_air_incident, absorb(RN_year t) cluster(RN_id )
unique t if e(sample)
replace b = _b[p_air_incident] if Years == `h'
replace se = _se[p_air_incident] if Years == `h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == `h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == `h'
}

// Regressions for pre-incident placebo outcomes
foreach h in 2 3 4 5 6 7 8 9 10 11 12 {
reghdfe p_air_noe_neg`h' p_air_incident, absorb(RN_year t) cluster(RN_id)
unique t if e(sample)
replace b = _b[p_air_incident] if Years == -`h'
replace se = _se[p_air_incident] if Years == -`h'
replace u = (_b[p_air_incident] + 1.96*_se[p_air_incident]) if Years == -`h'
replace d = (_b[p_air_incident] - 1.96*_se[p_air_incident]) if Years == -`h'

}


preserve
keep if Years != .
keep b u d se Years Zero
graph set window fontface "Times New Roman"
twoway(rarea u d Years, col(navy) fint(inten20) lwidth(0) lpattern(solid))(line b Years, lcolor(navy) lpattern(solid) lwidth(medium))(line Zero Years, lcolor(black)), xlabel(-12(1)12, nogrid labsize(vlarge)) ylabel(-.03(0.01)0.02, labsize(vlarge)) legend(off) ytitle("{&Delta} P(Notice of Enforcement)", size(vlarge)) xtitle("Month", size(vlarge)) graphregion(color(white)) plotregion(color(white)) xsize(8.6)

// Save figure
graph export "$figures/Figure_A4_Panel_C.pdf", replace

restore